#R version 4.2.1 (2022-06-23)
#Platform: x86_64-w64-mingw32/x64 (64-bit)
#Running under: Windows 10 x64 (build 19044)

library(dplyr) #version 1.0.9
library(readr) #version 2.1.2
library(stringr) #version 1.4.0
library(zoo) #version 1.8-10
library(sandwich) #version 3.0-2
library(lmtest) #version 0.9-40
library(prediction) #version 0.3.14
library(ggpubr) #version 0.4.0
library(estimatr) #version 1.0.0
library(margins) #version 0.3.26
library(broom) #version 1.0.0
library(scales) #version 1.2.0
library(stargazer) #version 5.2.3

`%notin%` <- function(x,y) !(x %in% y) 

#GWF Institutionalization Figure 3####
d1 <- read_csv("data/gwf_fa_pi.csv") 
d2 <- read_csv("data/gwf_irt.csv") %>% 
  rename(gwf_casename = case) %>%
  select(gwf_casename, transition, gwf_irt) 

gwf <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("gwf_casename", "year")) %>% 
  left_join(d2, by = c("gwf_casename")) %>% 
  group_by(gwf_casename) %>% 
  mutate_at(.vars = vars(gwf_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, gwf_casename, gwf_regimetype, v2xps_party, v2pssunpar, gwf_incumbent_pi, gwf_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, gwf_priorduration, coups, gwf_number, EFindex, theta_mean) %>%
  mutate(gwf_incumbent_pi = ifelse(is.na(gwf_incumbent_pi), v2xps_party, gwf_incumbent_pi)) %>% 
  rename(incumbent_pi = gwf_incumbent_pi) %>% 
  filter(transition == year, gwf_casename %notin% c("Bolivia 43-46", "Thailand 57-73", "Pakistan 58-71", "Haiti 88-90",
                                                    "Sudan 58-64", "Ecuador 63-66", "Venezuela 48-58")) %>% 
  ungroup() %>% 
  unique()

m1 <- lm_robust(gwf_irt ~ incumbent_pi, data = gwf, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Bivariate") 

m2 <- lm_robust(gwf_irt ~ incumbent_pi + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + gwf_priorduration + coups + gwf_number + EFindex + theta_mean, data = gwf, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Base") 

m3 <- lm_robust(gwf_irt ~ incumbent_pi + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + gwf_priorduration + coups + gwf_number + EFindex + theta_mean, data = gwf, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Inequality")

rbind(m1, m2, m3) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(conf.low = round(conf.low, 3), conf.high = round(conf.high, 3)) %>% 
  mutate(term = factor(term, levels = c("theta_mean", "EFindex", "gwf_number", "coups", "gwf_priorduration", "lag_milexp", "larea", "lgupop", "e_migdppcln", "e_regionpol", "e_total_resources_income_pc", "islam50", 
                                        "e_peaveduc", "ginibestestimateextrapolatednly", "LagRuralInequality", "incumbent_pi"))) %>% 
  mutate(Model = factor(Model, levels = c("Bivariate", "Base", "Inequality"))) %>% 
  mutate(term = recode(term, "theta_mean" = "Repression", "lag_milexp" = "Military Spending", "larea" = "Area", "lgupop" = "Urban", "e_migdppcln" = "GDP PC", "e_regionpol" = "Region", 
                       "e_total_resources_income_pc" = "Resources", "islam50" = "Islam", "e_peaveduc" = "Education", "ginibestestimateextrapolatednly" = "GINI", 
                       "LagRuralInequality" = "Rural", "incumbent_pi" = "Institutionalization", "gwf_priorduration" = "Duration", "coups" = "Coups", "gwf_number" = "Breakdowns", "EFindex" = "Ethnic Frac")) %>% 
  ggplot(., aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high, color = factor(Model), shape = Model)) +
  geom_pointrange(show.legend = F, position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept=0, linetype = "dashed") +
  geom_line(aes(y = estimate-100)) +
  geom_point(aes(y = estimate-100)) +
  theme_bw() +
  xlab("") +
  ylab("Bounded Democracy") +
  ylim(-2, 2) +
  scale_color_grey(name = "Model", start = 0.8, end = 0) +
  theme(legend.background = element_rect(fill = alpha(.1)), legend.position = "bottom") +
  coord_flip() 

ggsave(filename = "figures/fig3.jpg", dpi = 500, width = 5, height = 5) 

#GWF Strength models - Figure 4 #####
d1 <- read_csv("data/gwf_fa_strength.csv") 
d2 <- read_csv("data/gwf_irt.csv") %>% 
  rename(gwf_casename = case) %>%
  dplyr::select(gwf_casename, transition, gwf_irt) 

gwf <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("gwf_casename", "year")) %>% 
  left_join(d2, by = c("gwf_casename")) %>% 
  group_by(gwf_casename) %>% 
  mutate_at(.vars = vars(gwf_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

gwf$v2pssunpar <- scales::rescale(gwf$v2pssunpar, to = c(0, 1), from = range(gwf$v2pssunpar, na.rm = TRUE, finite = TRUE))

gwf %<>%
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, gwf_casename, gwf_regimetype, v2xps_party, v2pssunpar, gwf_incumbent_strength, 
                gwf_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, gwf_priorduration, coups, gwf_number, EFindex, theta_mean) %>%
  mutate(gwf_incumbent_strength = ifelse(is.na(gwf_incumbent_strength), v2pssunpar, gwf_incumbent_strength)) %>% 
  rename(incumbent_strength = gwf_incumbent_strength) %>% 
  #filter outliers
  filter(transition == year, gwf_casename %notin% c("Bolivia 43-46", "Thailand 57-73", "Pakistan 58-71", "Haiti 88-90",
                                                    "Sudan 58-64", "Ecuador 63-66", "Venezuela 48-58")) %>% 
  ungroup() %>% 
  unique() 

m1 <- lm_robust(gwf_irt ~ incumbent_strength, data = gwf, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Bivariate") 

m2 <- lm_robust(gwf_irt ~ incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + gwf_priorduration + coups + gwf_number + EFindex + theta_mean, data = gwf, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Base") 

m3 <- lm_robust(gwf_irt ~ incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + gwf_priorduration + coups + gwf_number + EFindex + theta_mean, data = gwf, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Inequality")

rbind(m1, m2, m3) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(conf.low = round(conf.low, 3), conf.high = round(conf.high, 3)) %>% 
  mutate(term = factor(term, levels = c("theta_mean", "EFindex", "gwf_number", "coups", "gwf_priorduration", "lag_milexp", "larea", "lgupop", "e_migdppcln", "e_regionpol", "e_total_resources_income_pc", "islam50", 
                                        "e_peaveduc", "ginibestestimateextrapolatednly", "LagRuralInequality", "incumbent_strength"))) %>% 
  mutate(Model = factor(Model, levels = c("Bivariate", "Base", "Inequality"))) %>% 
  mutate(term = recode(term, "theta_mean" = "Repression", "lag_milexp" = "Military Spending", "larea" = "Area", "lgupop" = "Urban", "e_migdppcln" = "GDP PC", "e_regionpol" = "Region", 
                       "e_total_resources_income_pc" = "Resources", "islam50" = "Islam", "e_peaveduc" = "Education", "ginibestestimateextrapolatednly" = "GINI", 
                       "LagRuralInequality" = "Rural", "incumbent_strength" = "Strength", "gwf_priorduration" = "Duration", "coups" = "Coups", "gwf_number" = "Breakdowns", "EFindex" = "Ethnic Frac")) %>% 
  ggplot(., aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high, color = factor(Model), shape = Model)) +
  geom_pointrange(show.legend = F, position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept=0, linetype = "dashed") +
  geom_line(aes(y = estimate-100)) +
  geom_point(aes(y = estimate-100)) +
  theme_bw() +
  xlab("") +
  ylab("Bounded Democracy") +
  ylim(-2, 2) +
  scale_color_grey(name = "Model", start = 0.8, end = 0) +
  theme(legend.background = element_rect(fill = alpha(.1)), legend.position = "bottom") +
  coord_flip() 

ggsave(filename = "figures/fig4.jpg", dpi = 500, width = 5, height = 5) 

#PAR Models Figure 5######
#institutionalization
d1 <- read_csv("data/svolik_fa_pi.csv") 
d2 <- read_csv("data/svolik_irt.csv") %>% 
  rename(svolik_case = case) %>%
  select(svolik_case, transition, svolik_irt)

svolik <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("svolik_case", "year")) %>% 
  left_join(d2, by = c("svolik_case")) %>% 
  group_by(svolik_case) %>% 
  mutate_at(.vars = vars(svolik_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, svolik_case, svolik_regime, v2xps_party, v2pssunpar, svolik_incumbent_pi, svolik_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea,
         svolik_priorduration, coups, svolik_number, EFindex, theta_mean) %>%
  mutate(svolik_incumbent_pi = ifelse(is.na(svolik_incumbent_pi), v2xps_party, svolik_incumbent_pi)) %>% 
  rename(incumbent_pi = svolik_incumbent_pi)  %>% 
  filter(transition == year, svolik_case %notin% c("Bolivia 1946-1946", "Thailand 1971-1972", "Pakistan 1958-1972", "Haiti 1986-1990", 
                                                   "Sudan 1959-1964", "Ecuador 1963-1965", "Haiti 1950-1956",
                                                   "Comoros 1999-2001", "Venezuela 1948-1949", "Burkina Faso (Upper Volta) 1966-1971", "Peru 1949-1956"))%>% 
  ungroup() %>% 
  unique()

m1 <- lm_robust(svolik_irt ~ incumbent_pi, data = svolik, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Bivariate") 

m2 <- lm_robust(svolik_irt ~ incumbent_pi + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = svolik, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Base") 

m3 <- lm_robust(svolik_irt ~ incumbent_pi + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = svolik, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Inequality")

rbind(m1, m2, m3) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(conf.low = round(conf.low, 3), conf.high = round(conf.high, 3)) %>% 
  mutate(term = factor(term, levels = c("theta_mean", "EFindex", "svolik_number", "coups", "svolik_priorduration", "lag_milexp", "larea", "lgupop", "e_migdppcln", "e_regionpol", "e_total_resources_income_pc", "islam50", 
                                        "e_peaveduc", "ginibestestimateextrapolatednly", "LagRuralInequality", "incumbent_pi"))) %>% 
  mutate(Model = factor(Model, levels = c("Bivariate", "Base", "Inequality"))) %>% 
  mutate(term = recode(term, "theta_mean" = "Repression", "lag_milexp" = "Military Spending", "larea" = "Area", "lgupop" = "Urban", "e_migdppcln" = "GDP PC", "e_regionpol" = "Region", 
                       "e_total_resources_income_pc" = "Resources", "islam50" = "Islam", "e_peaveduc" = "Education", "ginibestestimateextrapolatednly" = "GINI", 
                       "LagRuralInequality" = "Rural", "incumbent_pi" = "Institutionalization", "svolik_priorduration" = "Duration", "coups" = "Coups", "svolik_number" = "Breakdowns", "EFindex" = "Ethnic Frac")) %>% 
  ggplot(., aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high, color = factor(Model), shape = Model)) +
  geom_pointrange(show.legend = F, position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept=0, linetype = "dashed") +
  geom_line(aes(y = estimate-100)) +
  geom_point(aes(y = estimate-100)) +
  theme_bw() +
  xlab("Party Institutionalization") +
  ylab("Bounded Democracy") +
  ylim(-2, 2) +
  scale_color_grey(name = "Model", start = 0.8, end = 0) +
  theme(legend.background = element_rect(fill = alpha(.1)), legend.position = "bottom") +
  coord_flip() 

ggsave(filename = "figures/fig5_1.jpg", dpi = 500, width = 5, height = 5) 

#strength
d1 <- read_csv("data/svolik_fa_strength.csv") 
d2 <- read_csv("data/svolik_irt.csv") %>% 
  rename(svolik_case = case) %>%
  dplyr::select(svolik_case, transition, svolik_irt)

svolik <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("svolik_case", "year")) %>% 
  left_join(d2, by = c("svolik_case")) %>% 
  group_by(svolik_case) %>% 
  mutate_at(.vars = vars(svolik_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

svolik$v2pssunpar <- scales::rescale(svolik$v2pssunpar, to = c(0, 1), from = range(svolik$v2pssunpar, na.rm = TRUE, finite = TRUE))

svolik %<>%  
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar))%>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, svolik_case, svolik_regime, v2xps_party, v2pssunpar, svolik_incumbent_strength, 
                svolik_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea,
                svolik_priorduration, coups, svolik_number, EFindex, theta_mean) %>%
  mutate(svolik_incumbent_strength = ifelse(is.na(svolik_incumbent_strength), v2pssunpar, svolik_incumbent_strength)) %>% 
  rename(incumbent_strength = svolik_incumbent_strength) %>% 
  filter(transition == year, svolik_case %notin% c("Bolivia 1946-1946", "Thailand 1971-1972", "Pakistan 1958-1972", "Haiti 1986-1990", 
                                                   "Sudan 1959-1964", "Ecuador 1963-1965", "Haiti 1950-1956",
                                                   "Comoros 1999-2001", "Venezuela 1948-1949", "Burkina Faso (Upper Volta) 1966-1971", "Peru 1949-1956"))%>% 
  ungroup() %>% 
  unique()

m1 <- lm_robust(svolik_irt ~ incumbent_strength, data = svolik, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Bivariate") 

m2 <- lm_robust(svolik_irt ~ incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = svolik, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Base") 

m3 <- lm_robust(svolik_irt ~ incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = svolik, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Inequality")

rbind(m1, m2, m3) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(conf.low = round(conf.low, 3), conf.high = round(conf.high, 3)) %>% 
  mutate(term = factor(term, levels = c("theta_mean", "EFindex", "svolik_number", "coups", "svolik_priorduration", "lag_milexp", "larea", "lgupop", "e_migdppcln", "e_regionpol", "e_total_resources_income_pc", "islam50", 
                                        "e_peaveduc", "ginibestestimateextrapolatednly", "LagRuralInequality", "incumbent_strength"))) %>% 
  mutate(Model = factor(Model, levels = c("Bivariate", "Base", "Inequality"))) %>% 
  mutate(term = recode(term, "theta_mean" = "Repression", "lag_milexp" = "Military Spending", "larea" = "Area", "lgupop" = "Urban", "e_migdppcln" = "GDP PC", "e_regionpol" = "Region", 
                       "e_total_resources_income_pc" = "Resources", "islam50" = "Islam", "e_peaveduc" = "Education", "ginibestestimateextrapolatednly" = "GINI", 
                       "LagRuralInequality" = "Rural", "incumbent_strength" = "Strength", "svolik_priorduration" = "Duration", "coups" = "Coups", "svolik_number" = "Breakdowns", "EFindex" = "Ethnic Frac")) %>% 
  ggplot(., aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high, color = factor(Model), shape = Model)) +
  geom_pointrange(show.legend = F, position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept=0, linetype = "dashed") +
  geom_line(aes(y = estimate-100)) +
  geom_point(aes(y = estimate-100)) +
  theme_bw() +
  xlab("Party Strength") +
  ylab("Bounded Democracy") +
  ylim(-2, 2) +
  scale_color_grey(name = "Model", start = 0.8, end = 0) +
  theme(legend.background = element_rect(fill = alpha(.1)), legend.position = "bottom") +
  coord_flip() 

ggsave(filename = "figures/fig5_2.jpg", dpi = 500, width = 5, height = 5) 


#WTH Models Figure 6######
#institutionalization##
d1 <- read_csv("data/wth_fa_pi.csv") 
d2 <- read_csv("data/wth_irt.csv") %>% 
  rename(wth_case = case) %>%
  select(wth_case, transition, wth_irt)

wth <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("wth_case", "year")) %>% 
  left_join(d2, by = c("wth_case")) %>% 
  group_by(wth_case) %>% 
  mutate_at(.vars = vars(wth_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, wth_case, wth_regime, v2xps_party, v2pssunpar, wth_incumbent_pi, wth_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, 
         wth_priorduration, coups, wth_number, EFindex, theta_mean) %>%
  mutate(wth_incumbent_pi = ifelse(is.na(wth_incumbent_pi), v2xps_party, wth_incumbent_pi)) %>% 
  rename(incumbent_pi = wth_incumbent_pi)  %>% 
  filter(transition == year, wth_case %notin% c("Haiti 1986-1989", "Comoros 1999-2001")) %>% 
  ungroup() %>% 
  unique()

m1 <- lm_robust(wth_irt ~ incumbent_pi, data = wth, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Bivariate") 

m2 <- lm_robust(wth_irt ~ incumbent_pi + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + wth_priorduration + coups + wth_number + EFindex + theta_mean, data = wth, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Base") 

m3 <- lm_robust(wth_irt ~ incumbent_pi + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + wth_priorduration + coups + wth_number + EFindex + theta_mean, data = wth, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Inequality")

rbind(m1, m2, m3) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(conf.low = round(conf.low, 3), conf.high = round(conf.high, 3)) %>% 
  mutate(term = factor(term, levels = c("theta_mean", "EFindex", "wth_number", "coups", "wth_priorduration", "lag_milexp", "larea", "lgupop", "e_migdppcln", "e_regionpol", "e_total_resources_income_pc", "islam50", 
                                        "e_peaveduc", "ginibestestimateextrapolatednly", "LagRuralInequality", "incumbent_pi"))) %>% 
  mutate(Model = factor(Model, levels = c("Bivariate", "Base", "Inequality"))) %>% 
  mutate(term = recode(term, "theta_mean" = "Repression", "lag_milexp" = "Military Spending", "larea" = "Area", "lgupop" = "Urban", "e_migdppcln" = "GDP PC", "e_regionpol" = "Region", 
                       "e_total_resources_income_pc" = "Resources", "islam50" = "Islam", "e_peaveduc" = "Education", "ginibestestimateextrapolatednly" = "GINI", 
                       "LagRuralInequality" = "Rural", "incumbent_pi" = "Institutionalization", "wth_priorduration" = "Duration", "coups" = "Coups", "wth_number" = "Breakdowns", "EFindex" = "Ethnic Frac")) %>% 
  ggplot(., aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high, color = factor(Model), shape = Model)) +
  geom_pointrange(show.legend = F, position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept=0, linetype = "dashed") +
  geom_line(aes(y = estimate-100)) +
  geom_point(aes(y = estimate-100)) +
  theme_bw() +
  xlab("Party Institutionalization") +
  ylab("Bounded Democracy") +
  ylim(-2, 2) +
  scale_color_grey(name = "Model", start = 0.8, end = 0) +
  theme(legend.background = element_rect(fill = alpha(.1)), legend.position = "bottom") +
  coord_flip()  

ggsave(filename = "figures/fig6_1.jpg", dpi = 500, width = 5, height = 5) 

#strength models #
d1 <- read_csv("data/wth_fa_strength.csv") 
d2 <- read_csv("data/wth_irt.csv") %>% 
  rename(wth_case = case) %>%
  dplyr::select(wth_case, transition, wth_irt)

wth <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("wth_case", "year")) %>% 
  left_join(d2, by = c("wth_case")) %>% 
  group_by(wth_case) %>% 
  mutate_at(.vars = vars(wth_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

wth$v2pssunpar <- scales::rescale(wth$v2pssunpar, to = c(0, 1), from = range(wth$v2pssunpar, na.rm = TRUE, finite = TRUE))

wth %<>%
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, wth_case, wth_regime, v2xps_party, v2pssunpar, wth_incumbent_strength, 
                wth_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, 
                wth_priorduration, coups, wth_number, EFindex, theta_mean) %>%
  mutate(wth_incumbent_strength = ifelse(is.na(wth_incumbent_strength), v2pssunpar, wth_incumbent_strength)) %>% 
  rename(incumbent_strength = wth_incumbent_strength) %>% 
  filter(transition == year, wth_case %notin% c("Haiti 1986-1989", "Comoros 1999-2001")) %>% 
  ungroup() %>% 
  unique() 

m1 <- lm_robust(wth_irt ~ incumbent_strength, data = wth, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Bivariate") 

m2 <- lm_robust(wth_irt ~ incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + wth_priorduration + coups + wth_number + EFindex + theta_mean, data = wth, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Base") 

m3 <- lm_robust(wth_irt ~ incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + wth_priorduration + coups + wth_number + EFindex + theta_mean, data = wth, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Inequality")

rbind(m1, m2, m3) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(conf.low = round(conf.low, 3), conf.high = round(conf.high, 3)) %>% 
  mutate(term = factor(term, levels = c("theta_mean", "EFindex", "wth_number", "coups", "wth_priorduration", "lag_milexp", "larea", "lgupop", "e_migdppcln", "e_regionpol", "e_total_resources_income_pc", "islam50", 
                                        "e_peaveduc", "ginibestestimateextrapolatednly", "LagRuralInequality", "incumbent_strength"))) %>% 
  mutate(Model = factor(Model, levels = c("Bivariate", "Base", "Inequality"))) %>% 
  mutate(term = recode(term, "theta_mean" = "Repression", "lag_milexp" = "Military Spending", "larea" = "Area", "lgupop" = "Urban", "e_migdppcln" = "GDP PC", "e_regionpol" = "Region", 
                       "e_total_resources_income_pc" = "Resources", "islam50" = "Islam", "e_peaveduc" = "Education", "ginibestestimateextrapolatednly" = "GINI", 
                       "LagRuralInequality" = "Rural", "incumbent_strength" = "Strength", "wth_priorduration" = "Duration", "coups" = "Coups", "wth_number" = "Breakdowns", "EFindex" = "Ethnic Frac")) %>% 
  ggplot(., aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high, color = factor(Model), shape = Model)) +
  geom_pointrange(show.legend = F, position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept=0, linetype = "dashed") +
  geom_line(aes(y = estimate-100)) +
  geom_point(aes(y = estimate-100)) +
  theme_bw() +
  xlab("Party Strength") +
  ylab("Bounded Democracy") +
  ylim(-2, 2) +
  scale_color_grey(name = "Model", start = 0.8, end = 0) +
  theme(legend.background = element_rect(fill = alpha(.1)), legend.position = "bottom") +
  coord_flip()

ggsave(filename = "figures/fig6_2.jpg", dpi = 500, width = 5, height = 5) 

#DD Models Figure 7####
#institutionalization
d1 <- read_csv("data/dd_fa_pi.csv") 
d2 <- read_csv("data/dd_irt.csv") %>% 
  rename(dd_case = case) %>%
  select(dd_case, transition, dd_irt)

dd <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("dd_case", "year")) %>% 
  left_join(d2, by = c("dd_case")) %>% 
  group_by(dd_case) %>% 
  mutate_at(.vars = vars(dd_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, dd_case, dd_regime, v2xps_party, v2pssunpar, dd_incumbent_pi, dd_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, 
         dd_priorduration, coups, dd_number, EFindex, theta_mean) %>%
  mutate(dd_incumbent_pi = ifelse(is.na(dd_incumbent_pi), v2xps_party, dd_incumbent_pi)) %>% 
  rename(incumbent_pi = dd_incumbent_pi)  %>% 
  filter(transition == year, dd_case %notin% c("Bolivia 1946-1946", "Thailand 1946-1972", "Pakistan 1958-1970", "Haiti 1986-1989", 
                                               "Sudan 1958-1963", "Ecuador 1963-1965", "Haiti 1950-1955",
                                               "Comoros 1999-2003", "Venezuela 1948-1958")) %>% 
  ungroup() %>% 
  unique()

m1 <- lm_robust(dd_irt ~ incumbent_pi, data = dd, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Bivariate")

m2 <- lm_robust(dd_irt ~ incumbent_pi + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + dd_priorduration + coups + dd_number + EFindex + theta_mean, data = dd, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Base") 

m3 <- lm_robust(dd_irt ~ incumbent_pi + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + dd_priorduration + coups + dd_number + EFindex + theta_mean, data = dd, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Inequality")

rbind(m1, m2, m3) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(conf.low = round(conf.low, 3), conf.high = round(conf.high, 3)) %>% 
  mutate(term = factor(term, levels = c("theta_mean", "EFindex", "dd_number", "coups", "dd_priorduration", "lag_milexp", "larea", "lgupop", "e_migdppcln", "e_regionpol", "e_total_resources_income_pc", "islam50", 
                                        "e_peaveduc", "ginibestestimateextrapolatednly", "LagRuralInequality", "incumbent_pi"))) %>% 
  mutate(Model = factor(Model, levels = c("Bivariate", "Base", "Inequality"))) %>% 
  mutate(term = recode(term, "theta_mean" = "Repression", "lag_milexp" = "Military Spending", "larea" = "Area", "lgupop" = "Urban", "e_migdppcln" = "GDP PC", "e_regionpol" = "Region", 
                       "e_total_resources_income_pc" = "Resources", "islam50" = "Islam", "e_peaveduc" = "Education", "ginibestestimateextrapolatednly" = "GINI", 
                       "LagRuralInequality" = "Rural", "incumbent_pi" = "Institutionalization", "dd_priorduration" = "Duration", "coups" = "Coups", "dd_number" = "Breakdowns", "EFindex" = "Ethnic Frac")) %>% 
  ggplot(., aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high, color = factor(Model), shape = Model)) +
  geom_pointrange(show.legend = F, position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept=0, linetype = "dashed") +
  geom_line(aes(y = estimate-100)) +
  geom_point(aes(y = estimate-100)) +
  theme_bw() +
  xlab("Party Institutionalization") +
  ylab("Bounded Democracy") +
  ylim(-2, 2) +
  scale_color_grey(name = "Model", start = 0.8, end = 0) +
  theme(legend.background = element_rect(fill = alpha(.1)), legend.position = "bottom") +
  coord_flip() 

ggsave(filename = "figures/fig7_1.jpg", dpi = 500, width = 5, height = 5) 

#strength
d1 <- read_csv("data/dd_fa_strength.csv") 
d2 <- read_csv("data/dd_irt.csv") %>% 
  rename(dd_case = case) %>%
  dplyr::select(dd_case, transition, dd_irt)

dd <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("dd_case", "year")) %>% 
  left_join(d2, by = c("dd_case")) %>% 
  group_by(dd_case) %>% 
  mutate_at(.vars = vars(dd_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

dd$v2pssunpar <- scales::rescale(dd$v2pssunpar, to = c(0, 1), from = range(dd$v2pssunpar, na.rm = TRUE, finite = TRUE))

dd %<>%
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, dd_case, dd_regime, v2xps_party, v2pssunpar, dd_incumbent_strength, 
                dd_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, 
                dd_priorduration, coups, dd_number, EFindex, theta_mean) %>%
  mutate(dd_incumbent_strength = ifelse(is.na(dd_incumbent_strength), v2pssunpar, dd_incumbent_strength)) %>% 
  rename(incumbent_strength = dd_incumbent_strength) %>% 
  filter(transition == year, dd_case %notin% c("Bolivia 1946-1946", "Thailand 1946-1972", "Pakistan 1958-1970", "Haiti 1986-1989", 
                                               "Sudan 1958-1963", "Ecuador 1963-1965", "Haiti 1950-1955",
                                               "Comoros 1999-2003", "Venezuela 1948-1958")) %>% 
  ungroup() %>% 
  unique() 

m1 <- lm_robust(dd_irt ~ incumbent_strength, data = dd, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Bivariate")

m2 <- lm_robust(dd_irt ~ incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + dd_priorduration + coups + dd_number + EFindex + theta_mean, data = dd, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Base") 

m3 <- lm_robust(dd_irt ~ incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + dd_priorduration + coups + dd_number + EFindex + theta_mean, data = dd, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Inequality")

rbind(m1, m2, m3) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(conf.low = round(conf.low, 3), conf.high = round(conf.high, 3)) %>% 
  mutate(term = factor(term, levels = c("theta_mean", "EFindex", "dd_number", "coups", "dd_priorduration", "lag_milexp", "larea", "lgupop", "e_migdppcln", "e_regionpol", "e_total_resources_income_pc", "islam50", 
                                        "e_peaveduc", "ginibestestimateextrapolatednly", "LagRuralInequality", "incumbent_strength"))) %>% 
  mutate(Model = factor(Model, levels = c("Bivariate", "Base", "Inequality"))) %>% 
  mutate(term = recode(term, "theta_mean" = "Repression", "lag_milexp" = "Military Spending", "larea" = "Area", "lgupop" = "Urban", "e_migdppcln" = "GDP PC", "e_regionpol" = "Region", 
                       "e_total_resources_income_pc" = "Resources", "islam50" = "Islam", "e_peaveduc" = "Education", "ginibestestimateextrapolatednly" = "GINI", 
                       "LagRuralInequality" = "Rural", "incumbent_strength" = "Strength", "dd_priorduration" = "Duration", "coups" = "Coups", "dd_number" = "Breakdowns", "EFindex" = "Ethnic Frac")) %>% 
  ggplot(., aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high, color = factor(Model), shape = Model)) +
  geom_pointrange(show.legend = F, position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept=0, linetype = "dashed") +
  geom_line(aes(y = estimate-100)) +
  geom_point(aes(y = estimate-100)) +
  theme_bw() +
  xlab("Party Strength") +
  ylab("Bounded Democracy") +
  ylim(-2, 2) +
  scale_color_grey(name = "Model", start = 0.8, end = 0) +
  theme(legend.background = element_rect(fill = alpha(.1)), legend.position = "bottom") +
  coord_flip() 

ggsave(filename = "figures/fig7_2.jpg", dpi = 500, width = 5, height = 5) 
